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• Abstract 

■ The paper is concerned with the problem of a semi-infinite crack at the interface between two 

C I dissimilar elastic half-spaces, loaded by a general asymmetrical system of forces distributed along the 

. crack faces. On the basis of the weight function approach and the fundamental reciprocal identity 

(Betti formula), we formulate the elasticity problem in terms of singular integral equations relating 
, the applied loading and the resulting crack opening. Such formulation is fundamental in the theory 

of elasticity and extensively used to solve several problems in linear elastic fracture mechanics (for 
' (— I '■ instance various classic crack problems in homogeneous and heterogeneous media) . This formulation 

^ ^ is also crucial in important recent multiphysics applications, where the elastic problem is coupled 

I ■ with other concurrent physical phenomena. A paradigmatic example is hydraulic fracturing, where 

i^H ■ the elasticity equations are coupled with fluid dynamics. 
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1 Introduction 

CNJ , Three-dimensional problems of linear elasticity can be solved nowadays by a variety of techniques. Purely 

numerical methods, such as finite element and boundary element techniques, are very popular and 
extensively used, thanks to the computational power of modern computers. However, analytical methods, 
such as integral transforms, variational inequalities, and singular integral equations, are still of great 
, interest, for instance in deriving fundamental solutions for some basic geometries (to be employed in 

numerical methods), in the assertion of existence and uniqueness of the solution, and in bifurcation 
theory (Bigoni and Capuani, 2002; Bigoni and Dal Corso, 2008). 

In particular, the method of singular integral equations was first developed for two-dimensional prob- 
lems (Muskhelishvili, 1953) and later extended to three-dimensional problems by the rapidly developing 

■ theory of multi-dimensional singular integral operators (Kupradze et al., 1979; Mikhlin and Prossdorf, 
H , 1980). For crack problems, the general approach to derive the singular integral formulation is based 

■ on the Green's function method (Weaver, 1977; Budiansky and Rice, 1979, Linkov et al., 1997). This 
method consists in two steps. First, the stresses at a point x inside a body containing a discontinuity 
surface S are expressed in the form 

a^j{x)= [ K,jk{x,O^M^)d^, (1) 
Js 

where Kijk{x,^) = Cijst{d/dxt)T^i.{x,^) and T^f,{x,^) is the traction associated with the Green's func- 
tion Usp{x, ^), according to the equation r^"^(a;, 4) = n^CkipqUsp^q{x, ^). Here Cijst is the elasticity tensor, 
77,^ is the outward normal on the 4- side of the surface and the displacement discontinuity is denoted by 
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Aufc = Uf, — u^. The second step consists in taking the Umit as the point x approaches a point Xq on 
the surface S, to obtain the tractions on the surface 

p+(a;o)= lim f n+{x)K,jkix,^)Aukmt (2) 

The described procedure works when the tractions on the discontinuity surface 5* are symmetrical. 
Moreover, the Hmit in ^ is not straightforward and due care should be taken to ensure the convergence 
of the integrals. 

In this paper, we consider the three-dimensional problem of a semi-infinite crack at the interface 
between two dissimilar elastic half-spaces and apply a method which avoids the use of the Green's 
function and the tedious limiting procedure. We also do not use the assumption of symmetrical load. 

Our results are based on integral transforms and two fundamental notions of linear elasticity: the 
Betti reciprocal theorem, on the one hand, and the weight function approach, on the other. The Betti 
identity has been extensively used in the perturbation analysis of two and three-dimensional cracks 
(Willis and Movchan, 1995; Bercial-Velez et al., 2005). The concept of weight function was introduced in 
linear fracture mechanics by Bueckner (1970) for two-dimensional problem and later extended to general 
three-dimensional problems by Willis and Movchan (1995), Piccolroaz et al. (2007). 

We consider a semi-infinite crack at the interface between two dissimilar elastic half-spaces Rj_ = 
{x = (xi,X2,X3) £ R"^ : ±X2 > 0}, see Fig. [TJ The two elastic half-spaces are assumed to be isotropic 
with shear modulus and Poisson's ratio denoted by fi± and v±, respectively. The crack lies on the half- 
plane ~ {x — {xi,X2,X3) £ M.^ : xi < 0,X2 = 0}. The crack faces are loaded by a system of, not 
necessarily symmetrical, distributed forces p^. . 



















x^ 
V_, |I_ 



Figure 1: Semi-infinite interfacial crack loaded by asymmetrical forces on the crack faces. 

It is convenient to introduce the symmetrical and skew-symmetrical parts of the loading as follows 

(p> = ^(p+ + P-), bl=p+-p-, (3) 

where we used standard notations to denote the average, (/), and the jump, [/]], of a function / across 
the plane containing the crack, X2~^; 

{f)ixi,X3) - i[/(xi,0+,x-3) + /(a;i,0-,X3)], [/Ka:i, X3) = /(xi, 0+, zg) - /(xi, 0", xg). (4) 

In Section [21 we introduce the fundamental reciprocal identity and the concept of weight function, as 
special singular solution to the homogeneous crack problem. These preliminary results are then used 
in the sequel for the analysis of both two and three-dimensional elastic problems. In particular, the 
two-dimensional problem is analysed in Section [3] (Mode III) and Section H] (Mode I and II), while the 
general three-dimensional probelm is analysed in Section [5) 
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2 Preliminary Results. The Betti formula and weight functions 



In this section, we extend the Betti formula to the case of general asymmetrical loading applied at 
the crack surfaces. The Betti formula is used in order to relate the physical solution to the weight 
function, which is a special singular solution to the homogeneous problem (traction-free crack faces) 
(Willis and Movchan, 1995; Piccolroaz et al.. 2007). For the sake of generality, we allow for traction and 
displacement jumps across the interface, so that results are general and applicable also to the case of 
imperfect interfaces. 

In the absence of body forces, the Betti formula takes the form 

{cr(i)n.tt(2) -cr(2)n.„(i)}ds = o, (5) 

an 

where dQ is any surface enclosing a region Q within which both displacement fields it^^' and u^^^ satisfy 
the equations of equilibrium, with corresponding stress states cr'^^ and cr^^\ and n denotes the outward 
normal to dil. 

Applying the Betti formula to a hemispherical domain in the upper half-space R'^, whose plane 
boundary is X2 — 0^ and whose radius R will be allowed to tend to infinity, we obtain, in the limit 
R—T'oo, 

J {(Ti'\xi,0+, xs) ■ u^^'>{xi,0+, X3) ~ cr^')(xi, 0+, xa) • ^(^^(xi, 0+, X3)}rfa:irfx3 = 0, (6) 
(2:2=0+) 

provided that the fields u^^^ and it^^^ decay suitably fast at infinity. The notation (T2 is used to denote 
the traction vector acting on the plane a;2 = 0: ct2~ ce2- 

We can assume that u^^^ represents the physical field associated with the crack loaded at its surface, 
whereas u^"^^ represents a non-trivial solution of the homogeneous problem, the so-called weight function, 
defined as follows 

U^'^^{xx,X2,X-i) = RU{-Xi,X2,-X3), (7) 

where i? is a rotation matrix 

'-10 " 

i2= 1 . (8) 

-1 _ 

Note that the transformation ([7|) corresponds to introducing a change of coordinates in the solution it*^^' , 
namely a rotation about the a;2-axis through an angle tt. It is straightforward to verify that the weight 
function U satisfies the equations of equilibrium, but in a different domain where the crack is placed 
along the semi-plane X2 = 0, xi > 0. The notation S will be used for components of stress corresponding 
to the displacement field U, 

CT^'^^Xi,X2,X3) = RT,{-Xi,X2, ''X3)R. (9) 

Replacing u^'^^ {xi, X2, X3) with u^^\xi — x'i,X2,X3 — x'3), which corresponds to a shift within the plane 
{xi,X3), we obtain 

J {RU{x[ -si, 0^,X3 -xa) • (T2{xi,0~^ , X3) - R1^2ix'i - xijO'^ ,x'3- X3) ■u{xi,0^ ,X3)}dxidx3 = 0. 
(2:2=0+) 

(10) 

A similar equation can be derived by applying the Betti formula to a hemispherical domain in the lower 
half-space R'L , 

J {RU{x[ 0^, 0:3 -X3) ■ cr2{xi,0~ , X3) - RE2{x'i - xi,0^ , X3- X3) ■ u{xi,0~ , X3)}dxidx3 = 0. 
(2:2=0-) 

(11) 
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Subtracting ([TI]) from PU)) . we obtain 

/ {RlUjix\ - .xi,X3 - 0:3) • (cr2)(xi,.T3) + R{U){x\ - xi,x'3- X3) ■ [0-2] (xi, 2:3)- 
(.io) (12) 

-i?(S2)(xi - a;i,a;3 - X3) • a;3)} dzidzs = 0. 

Let us introduce the notations 

/( + )(xi,X3) = /(rEl,X3)7J(xi), /(-)(xi,X3) = /(xi,X3)i?(-Xi), (13) 
where H denotes the Heaviside function, so that 

fixi.Xs) ^ f<^+\xi,X3) + f^-\xi,X3). 

The reciprocity identity (fT^ becomes 

/ ■! filL^Ka:;'! -xi,a:;3 ~ X3) ■ ((T2)'+^(a;i, 2:3) -ii(S2)(xi -a;i,a;3 -0:3) • luf^\xi,X3)'>dxidx3 = 
J —00 

- / I ii[J7|(a::'i - xi,x'r^ ~ X3) ■ {(T2)'-~Hxi, X3) + R{U){x[ - xi,x'^ - X3) ■ l(T2i'-~\xi, X3) \ dxidxs 

- f {R{U){x[ ~XI,X'3- X3) ■ I(T2l(+)(xi,X3) - fi(S2>(a;i - Xi,.T^ - X3) ■ lui^+\xi,X3)] dxidX3. 

(14) 

The identity (|14p can be written in an equivalent form using the convolution with respect to both xi 
and X3 variables, denoted by the symbol ®, 

RlUl®{(T2)^+^ -R{'Z2)®luf'^ = -RlUl®{(T2)^-^-R{U)®l(T2l^~^-R{U)®l<T2l^+^+R{'S2)®lui<-+\ 

(15) 

Note that the superscripts '^^^ and ' denote functions whose support is restricted to the positive 
and negative semi-axes, respectively. Thus, (0-2) is the traction along the interface, ahead of the crack 
tip, whereas M(-) is the crack opening (displacement discontinuity across the crack faces). 

In the sequel, we will assume perfect contact conditions at the interface, so that traction and dis- 
placement are continuous for xi > 0: |(T2]]'^'' — ImJ'^-' — 0. Finally, the symmetric and skew-symmetric 
loads applied on the crack faces will be denoted by ((T2)'~^ = {p) and |(T2]]'~' = [pj, respectively. 

The terms |?7]] and (U) are known as the symmetric and skew-symmetric weight functions. The 
notion of skew-symmetric weight function was introduced in Piccolroaz et al. (2009) for problems of 
cracks loaded by asymmetrical forces. The term (S2) stands for the traction along the plane X2 ~ 0, 
corresponding to the singular solution U. 

3 Mode III 

In the case of antiplane deformation, the Betti formula psp relating the physical field U3, (T23 with the 
weight function U3, S23 reduces to the scalar equation 

IC/3I * (a23>(+) - (S23) * lu3f-^ = -IC/3I * (P3) - (C/3) * [P3l, (16) 

where the symbol * denotes the convolution with respect to the variable xi. 

Let us introduce the Fourier transform with respect to the variable xi as follows 

/oo 
/(xi)e"^-^dxi. (17) 
-00 

Applying the Fourier transform to the identity (|16|) we get 

1^31(^23)+ - (S23)Iu3r - -I^^31(P3) - (^73)1^31, (18) 
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where the superscripts and denote functions analytic in the upper and lower half-planes, respectively. 
We can now multiply both sides of ([T5)) by [C/s]]"^ to obtain 



where the factors in front of unknown functions are given by 



(19) 



(20) 



They can be computed from the relationships, which hold for the symmetric and skew- symmetric weight 
functions (Piccolroaz et al., 2009) 



b + e 



Thus, we can easily obtain 



A = 



V 



(S23), (f/3) 



B = -- 



2' h + e 

where h, e and are the following bimaterial constants, 



b = 



e = H , r] = ■ . 



(21) 



(22) 



(23) 



If we apply the inverse Fourier transform to (|19p . we obtain for the two opposite cases xi < Q and 
xi > the following relationships: 



-1 

2;i<0 



(24) 



(25) 



Note that the term ((723)^ in cancels from because it is a function, while the terms (pa) 
and [pa] cancel from (pS)) because they are "— " functions. 

Finally, we find the inverse Fourier transform of the function |/?||{t3|^. To this purpose we write 
1/^111^31^ ~ sign(/?) • /3|{{t3]~ and observe that sign(/3) is the Fourier transform of — i/(7ra;i), whereas 
/SJus]]^ is the Fourier transform of id\u3\~ / dxi. Therefore we get 



"3j 



1 , ^M' 

TTXl dxi 



(26) 



Note that, with slight abuse of notation, we write the convolution of two functions ,f{xi) and g{xi) in 
the form 

/•OO 

fi^i-09{0d^^f{xi)*g{xi), (27) 
so that the right-hand side in ([25)1 has to be interpreted as follows 

1 dlu3^- 



1 ^dju.l 

TTXl dxi 



(28) 



To simplify notations, we introduce the singular operator S and the orthogonal projectors V± {V++V- 
T) acting on the real axis: 



ip ^Sip^ * ip[;xi) = - / — 

TTXl TT Xl 



f{xi) ±xi > 0, 
otherwise. 



(29) 
(30) 
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Note that the operator S is the classic singular operator of the Cauchy type; it transforms any function 
if, satisfying the Holder condition, into a new function Stp which also satisfies the Holder condition 
(Muskhelishvili, 1946). The properties of the operator S in other functional spaces have been discussed 
in Prossdorf (1974). 

The integral identities (j24|) and ((25l) for Mode III deformation become 



W + fM = -^5(^)«^, ..<o, (31) 

2 + e oxi 

(a.3>W^--^5(^)«:l, ,,>o, (32) 
+ e oxi 

where S'^"'^ = VSV- is a singular operator, whereas iS^^-* = V+SV- is a compact operator (see Gakhov 
and Cherski 1978; Krein, 1958; Gohberg and Krein, 1958). Despite the fact that these operators look 
similar at a first glance, they are essentially different. Indeed, S'^'^^ : F(M_) F(R^), while iS*^^-* : 
F(M_) — > F(R+), where F(K.±) is some functional space of functions defined on M±. 
To make this point clear, let us write eqs. ([3T|) and ([5^ in extended form as 

(P3 + xbsl = 7I-— T / 7 — — d£„ xi < 0, (33) 

1 r ^a^t!*. (34, 

7r(6 + e) - ^ 

The integral in p3p is a Cauchy singular integral with a moving singularity, whereas the integral in (|34p 
is an integral with a fixed point singularity. 

Note that, in the particular case of a homogeneous body, we have 77 = and b + e = '^/fJ-, so that the 
integral identities (PTjl and ((5^ reduce to 

(R3)=-f5(^)^I|^, .KO, (35) 
2 axi 

(a.3>W=-^5(^)^^, XI >0. (36) 

The two equations ([55]) and form the system of integral identities for the Mode III deformation. 
The first equation p3p provides the integral relationship between the applied loading (^3), Ipsl and 
the resulting crack opening [[w3]](~-'. This is a singular integral equation and it is, generally speaking, 
invertible. However, the inverse operator depends on the properties of the solution. The second equation 
can be considered as an additional equation which allows to define the proper behaviour of the 
solution |u3|'-^'' and also, after the first equation has been inverted, to evaluate the traction ahead of the 
crack tip (cr23)*-~'^^- The reason for this is that the operator in the right-hand side of ([34]) is a compact 
one and thus it is not invertible. 

More details on the theory of singular integral equations can be found in Muskhelishvili (1946). The 
inversion of the singular operator 5'-''' , in some specific cases of interest for fracture mechanics and, in 
particular, for hydraulic fracturing, is discussed in Appendix 1X1 

4 Mode I and II 

In the case of plane strain deformation, the Betti identity (fT5)) relating the physical solution u = [ui,U2]'^ , 
<T2 = [c2i, CT22]"^ with the weight function U, S2 is given by 

nm * ('^2>(+' - i?(S2) * M^"^ = -mm * (p) - r{u) * m. m 
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where (p) = [{pi), (^2)]"^, M = [[pi], [P2l]^ are the symmetric and skew-symmetric parts of the loading. 
Again we assume ideal contact conditions, so that traction and displacement are continuous across the 
interface: [[cr2]]^+-' = [m]^^-' = 0. Here and in the sequel of this section, we use the following matrices: 



R = 



" -1 ■ 




■ 1 


■ 


, E = 


■ 1 ■ 


1 


, / = 





1 


-1 



(38) 



Note that the symmetric and skew-symmetric weight functions, |{[7| and (U), and the corresponding 
traction (S2) are represented by 2x2 matrices. In fact, in the case of an elastic bimaterial plane, there 
are two linearly independent weight functions, U-' = [?7j', [/j]"^: ^2 ^ [^21' ^22]'^; j — Ij^, and it is 
possible to construct the weight function tensors by ordering the components of each weight function in 
columns of 2x2 matrices (Piccolroaz et al., 2009) 



U 



^21 ^21 
^22 ^22 



Applying the Fourier transform to the equation (|37p . we obtain 

[«7f i2(^2)+ - {^2fRM- = -lufR{p) - {ufRM- 

Multiplying both sides by we get 

where A and B are the following matrices 



A = R-'lur^{UfR, B = R-'lur^{t2fR, 



(39) 



(40) 



(41) 



(42) 



which can be computed using results for the symmetric and skew-symmetric weight functions obtained 
in Antipov (1999) and Piccolroaz et al. (2009). Namely 



1 



lUl^-—[bI-zdsigniP)E] (S2 



{if) = ~— [aI - i-/sign{l3)E]{t2), 



(43) 



where b is defined in (j23p . d/b and a are the so-called Dundurs parameters, while 7 is another bimaterial 
constant: 



2m+ 



1 - 2v^ 



2/_i_ ' 

M-(l-2;/+) + /,t+(l-2;/_) 
2^_(1 - i/+) + 2/i+(l - v^y 



(44) 
(45) 



As a result, we obtain 



A = 



B 



[bl + ids\gTL{fi)E\ 



52-^2 

Inverting the Fourier transform for the two cases xi < and xi > 0, we get 



Alp 
Alp 



-1 

2-1 <o 



-1 

a:i>0 



Blu\- 

Bm 



(46) 
(47) 

(48) 
(49) 



Similarly to the previous section, the term (cr)+ in ([3D]) cancels from (pS]) because it is a function, 
whereas the term (p) cancels from (|49| because it is a "— " function. 
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To proceed further, we need to invert functions of the form sign(/3)|pi], and This is 

done in a way similar to the previous section and we finahy obtain the integral identities for plane strain 
deformation as follows 



dxi 
dxi 



xi < 0, 



xi > 0, 



(50) 



(51) 



where B^'*^ : F{1 



i^(R_), and A^''\B^''^ : F{R^) F(M+) are the following matrix operators 



2(62 - fi2) 



{ba- d-f)I + {da-b-f)ES'^'''>] , 5^")=- 



1 



b{da — 67) 



62 -d2 



62-^2 



blS'^"^ - dE 



(52) 



(53) 



2(62 _ ^2) 

Equations ([50]) and ([?T|) . together with the definition of operators ([5^ and form the system of 

integral identities for Mode I/II deformation. In particular, (|50p is a system of two singular integral 
equations and this system is coupled, in general. 

Note that, in the case where the Dundurs parameter d vanishes, the operators can be written as 



2 2 6 
2 6 



This means that the system ([50]) decouples in the case d = and reduces to 

xi < 0. 



6 0x1 



2 2 



(54) 
(55) 

(56) 



In the particular case of a homogeneous material, we additionally have a = 0, 6 = 2(1 — v)/^ and 
7 = (1 — 2z^)/[2(l — v)]^ so that the matrix operators simplify further to the form: 



^ " 4(1 -t.) '^ ' ^ ~ 2(l-t.) 



and the singular integral equation (|56p simplifies to 



2 



(-) 



{l-i^){p)-^-^ES^^\pl XI <0. 



(57) 



(58) 



In the case of symmetrical load {\p\ = 0) this formula can be found for example in Rice (1968). It differs 
from the Mode III case ([55]) only for the coefficient in front of (p). If in addition = 0, then the two 
formulas fully coincide, altough one is scalar and the other is vector one. 

The solution of either equation ([55)1 or ([55]) , allowing for the calculation of the crack opening related 
to a given load applied on the crack faces, requires the inversion of the operator S^'^\ which is discussed 
in details in Appendix [A] To invert the operator ((50)) in the general case, one needs to use results from 
Vekua (1970) where the systems of singular integral equations are discussed and analysed thoroughly. 

However, we propose here a simple way to decouple the system (j50p by acting on it by the following 
operator: 

n = 6/5(") + dE. (59) 
As a result, equation ([50)1 decouples as follows 



1 - P^)^^^ = \ [balS^'^ + E [67 +p(da - b^)K]^ lp\ + (blS^'^ + dE^ (p), 



(60) 
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where we have introduced a new bimaterial constant (0 < p < 1) 

In deriving (j60p we used results from Muskhehshvih (1946) to prove that 

(s'-'^y = -I + JC, (62) 
where X is the identity operator and K, : F(M_) F(R^) is a compact operator given by: 

1 f log|C/a;i| 
= ^ / ^'^(C)t^C- (63) 

It is easy to see that in the case of a homogeneous body (d = 0, a = and thus p = 1) the equation (pO|) 
is equivalent to (j58p . In the general case, the operator 1 — pK, can be analytically inverted again using 
results from Muskhehshvih (1946). Alternatively, the Fredholm equation ([BUI can be easily inverted by 
numerical techniques. 



5 General 3D case 



For the case of an elastic bimaterial space, there are three linearly independent weight functions, from 
which we can construct the weight function tensors by ordering the three components of each weight 
function in columns of 3x3 matrices (Piccolroaz et al., 2007; Piccolroaz et al., 2011) 



(64) 





\ ^1 


ul 


ul 1 




r yi 


^21 


y3 
^21 


u = 










^22 


y2 
^22 


y3 
^22 




ul 


ui 


ui 




^23 


^23 


^23 



Assuming again perfect interface (so that [cr- 



= M 



|(+) = 



= 0) , the Betti identity, relating the physical 



solution u = [ui, M2, 1*3]"^, cr2 = [cr2i,cr22, 7 0'23]^ with the weight function [/, S2, writes 

Riui ® (tT2)^+) - i2(S2) ® M^"^ = -RWl ® (p) - R{U) ® IpI 



(65) 



where the symbol ® denotes convolution with respect to both variables xi and X3 and R is the matrix 
defined in Sec. [T] 

Let us introduce the Fourier transform with respect to the xi and 2:3 coordinates as follows 



/(/3,a;2,A) 



/(rEi,X2,X3)e^^"i+'^^^^d:Eidx3. 



— 00 ^ — 00 



Taking the Fourier transform of (|65p we obtain 

VufR{^2)+-{^2fRM- = - 

Multiplying both sides of ^ by i?"^|C7]l"'^ we get 

(p)-+A[pr = -(^2 

where 



lUf R{p) - {Uy RM 



A = R-^IU\-^ iUy R, B = R-^IUY 



R. 



(66) 

(67) 

(68) 
(69) 



The two matrices A and B can be obtained from the following relationships between symmetric and 
skew-symmetric weight functions (Antipov, 1999; Piccolroaz et al., 2011) 



m = iG(S2) 



{U) = iF(E2), 



(70) 
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where p = + A^, 



G(/3, A) = - — [bp^I + eE2 + idpE^) , A) = - — (6 V/ + je-^ + z67p£;: 





" A2 





-p\ ' 




" 


-/3 


" 


^^2(/3,A) = 











, £^3(/3,A) = 


/? 





A 




-/3A 












-A 






I is the 3x3 identity matrix, E2E3 = and e, / are the following bimaterial parameters: 

1^+ V- 



e = — + — 



/ = 



(71) 
(72) 

(73) 



As a result, we obtain 



A(/3,A) 



2(62-(i2)(5 + g)^2 

{6(6a - d7)(6 + e)p2j + [6^7(6 + e) - 5a(6e + + (6^ - d^)]\ E2 + ^6(67 - da){b + e)p£;3}, 



(74) 



(75) 



By inversion of the Fourier transform for the two cases xi < Q and xi > 0, we obtain the two formulae 



-p-1 -p-l 
xi>0-' 3:3 



'aIp 



_ 1 1 



_ 1 1 



BM- 

BVa] 



(76) 
(77) 



We need now to perform the double Fourier inversion. The details are given in Appendix [E] The final 
result is given in terms of the following operators 



1 , , 1 



00 poo 



^(^1,6) 



1 ip{xi,X3) 



1 



00 nOO 



-00 V(2^1-6)2 + (x3-6)' 
{Xj -CjMCl,^3) 



-00 {xi-^i)^ + (2:3-6)^ 



n{xl+xl) 

The integral identities for 3D elasticity are given by 

{p)+A^^^M^B^^^luf-\ XI <0, 

(rT2) ( + ) + H = B^-^) M , XI > 0, 



d^id^3, i = l,3. 



(78) 

(79) 
(80) 

(81) 
(82) 



where 



2(62-d2)(6 + e) 
l^b{ba - dj){b + e)I+^ [bdj{b + e) ~ ba{be + d^) + {h^ ~ d^)f] £ 



2 



(83) 



1, 



-b{bj - da){b + e)Q^'^ E: 



d d 
dxi ' dx^ 
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b{b 



2(62 



d^){b + e) 



dxl 



ibe + d^)Q^'^E2 



d d 
dxi ' dx3 



2d{b + e)E3 



d _d_ 

dxi ' dx3 



(84) 



-[bdj{b 



1 

2{b^ - d^){b + e)^ 

) - ba{be + d^) 



(6^-.V]5r-i5(5.-.«)(5 + e)Q(^)^3(^,^ 



(85) 



1 



2{b^ - d^){b + e) 



h{h 



and the matrix operator E^^''^^ is defined as 



.(s,c) 







(s,c) 



- 9 - 




T 















d 






dxi . 







' dx3 



-Q 







(s,c) 



{s,c) 



9X3 



9 



Q\ 







9 



dxi " 9xi . 

In the case of homogeneous body the matrix operators A!''^\ S'*', S'-^-' simplify to 



8(1-:.) 



4(1 - z.) 



92 ^ 92 



9x2 





92 



dxidxs 






9 
9xi 





9 





9 


9xi 


9X3 





9 

9X3 







92 


92 




dxl 


dxj 









(1- 



92 



dxidxs 




9x2 



9x2 



(86) 



(87) 



(89) 



We note that the singular integral formulation (j8ip for the case of a homogeneous body and symmet- 
rical forces (Jp] = 0) is equivalent to the formulation given by Weaver (1977) and Budiansky and Rice 
(1979). To show this it is sufficient to use integration by parts, taking into account the behaviour of the 
displacements jump at infinity. 



6 Conclusions 

The singular integral formulation for a three-dimensional semi-infinite interfacial crack has been derived 
in explicit form by means of fundamental properties of linear elasticity and integral transforms. This 
formulation for interfacial cracks seems to be unknown in the literature. The method of derivation avoids 
the use of Green's function and the tedious limiting process involved in the general procedure. The 
identities are important in applications especially in case of multiphysics problems where the elasticity 
equations are coupled with other concurrent phenomena, such as in hydraulic fracturing. 
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A Inversion of the 2D singular operator S^'^^ 

The inversion formula for the singular integral equation 

V'(xi) = = - r ^^rf?, < 0, (A.l) 

strongly depends on the properties of the known function ip{xi). Indeed, let us first assume that the 
function tpixi) has compact support —a < xi < —6, where a and b are positive constants, and belongs 
to a Holder class, as it usually takes place in the classical 2D elasticity (Muskhelishvili, 1946). Then the 
inversion formula can be found in Muskhelishvili (1946) (see also Rice, 1968) 

^(xi) = (5(-^))"'v(xi)--i r ./^^^d^, xi<0. (A.2) 

V y 1 TT \ XlXl-^ 

Note that under such assumptions 

ip{xi) ^ a-i-s-0", ip{xi) ^ ^"Tc,/.. , xi^-oo. (A. 3) 

V-xi (-xi)-^/^ 

Here the constants Kq and K^o {Ko is the so-called stress intensity factor) are determined by the formulae 

K^ = --f ^d^, K^^- f mV~idt (A.4) 



Let us now consider the case where the function ip{xi) extends over the whole negative semi-axis, 
having the following behaviour at zero and infinity: 

~ , xi^O", ip{xi) ^ , xi^-CG. (A.5) 

(-Xi)"0 (-Xi)"°° 

If we assume that < 1/2 and aoo > 3/2 then the inversion formula (jA.2[) is still valid and leads to 
the asymptotics (jA.Sp . In the case ao < 1/2 and 1/2 < aoo < 3/2 the inversion formula (|A.2p remains 
true, the behaviour of the function 'ipi^i) near zero is the same as in (jA.3P i . while the behaviour at 



infinity (|A.3P 9 changes to: 

. ^ -7r'0ooSin(7rQ;oo) ^ . . 

(-a;i)"~ 

However, there are situations (such as some problems arising in hydraulic fracturing, see Garagash 
and Detournay, 2000) where the behaviour of the function "0(2^1) infinity is worst, so that aoo hi 
(jA.5p is in the range < aoo < 1/2. In such cases, the classic inversion formula (jA.2p is not any longer 
applicable and another inversion formula should be used instead. In Garagash and Detournay (2000) the 
formula 
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was derived by easy and elegant arguments leaving, however, the constant Co unknown. 

Here we present an inversion formula accurately derived under slightly more restricted conditions on 
the behaviour of the function ipi^i) B^t infinity. Namely we assume that 



1p{xi) = V'oo(2:i) + O 



xi — oo. 



where 



N 



U) 



(A., 



(A.9) 



and < a"i^ < 1/2 and j3 > 1/2. The alternative formula reads 



-1 ^ 
^(xi) = i^{x,) - ^tan(7ra(i)) 



(j) 



{-x,y^ 



Xl- ^ 



d^, xi < 0. (A.IO) 



Note that this formula can be rewritten in the equivalent form (|A.7p . where the constant Co can be 
determined by 



Co 



1 r ^(g)-V'oo(e) 



(A.ll) 



B Reduction of the identities from 3D to 2D case 



The integral identities for 3D elasticity can be reduced to the 2D case by assuming that all the mechanical 
fields are independent of the X3 coordinate. We start by showing that if the fields u, (T2, p are independent 
of X3 then the convolution with respect to Xi and X3, denoted by ®, reduces to the convolution with 
respect to the xi coordinate only, denoted by *. 

In particular, we have the following formulae (proved in Appendix[Dl Theorem lD.ip . If 93 is a function 
of xi only, then 

1 d'^ip ^ I dip 

dx\ 



\/xl+: 



xi dxi 



ox'l 0x1 



dxi xi dxi 



-2Sip. 



(B.l) 
(B.2) 



We also make use of the following identities (proved in Appendix [Dl Theorems ID. 21 and lD.3[) : 

xi dip X3 



x{ + X3 dxi x{ + xl dxs 
X3 dip x\ 



2 oT" ~ ^"""^ 



dx\ dx3 



dip 
dx-i 



dp 
dxi 

dip 



If iy9 is a function of x'l only then the identities (|B.3|) and (jB.4p yield 

xi dip 



„ y. dip 
dxi 0x3 



X ~\~ X 2 



2:3 



dxi 

dip 
dxi 



lixip 



® 



dpi 
dx\ 



(B.3) 
(B.4) 

(B.5) 
(B.6) 



Using the formulae (|B.1|) . (|B.2p . (jB.5|) and (jB.6|) . the integral identity (|5T|) for 3D elasticity reduces 



to 



2(62 „ ^2) 



(6a - d7) (doi - 67)5''*) 
-{da-h-i)S^'^ (ba-d-f) 








' 





" 


H- 























1 



(B.7) 
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1 



dxi 



dxi 







bSi^^^ 
oxi oxi 


















S^"^ 



d 
dxi 



(B.8) 



This formula is fully consistent with the formulation for 2D elasticity given in Sees. [3] and 01 In a similar 
way, also the formula ((82)) can be shown to reduce to the 2D formulation. 



C Illustrative example in 2D case: point forces applied at the 
crack faces 

In this section, we illustrate the use of the singular integral identities in the 2D case, when the loading 
is given in terms of point forces applied at the crack faces. We analyse first the antiplane problem and 
then the vector problem. 

Mode III. Symmetrical point forces. Assume that the loading is given as two symmetrical point 
forces applied on the crack faces, at a distance a from the crack tip, and directed along the zs-axis. 

(p3)(a:i) = -F5(xi+a), ipz\{x^) = Q, (C.l) 

where 5 is the Dirac delta function. 

The singular integral equation relating the applied loading with the resulting crack opening is given 
by (j3ip and, by means of the inversion formula \K.2\ we obtain 



dxi TT J-oo \ Xl Xl ~ ^ TT \ Xl Xl + a 



By integration of (jC.2p and using the conditions that the crack opening vanishes at zero and at infinity, 
we get 

'2F „ ^ , rix^ 

(b + e) arctanh W , — a<a;i<0, 

M(xi)=<( " (C.3) 
(o + e) arctanh . / , xi < —a. 

TT \ Xl 

From p2p we can now obtain the corresponding tractions ahead of the crack tip, namely 



(^23>'+H-^i) = -^J—-, / 7^^4^d^ = -\ — ■ (C.4) 

7r(o + e) j„oo — 4 ot, Tr\jxiXi+a 

Note that the stress intensity factor in this case is 

Km = lim V^^M^+H^i) = \ —F- (C.5) 

Mode III. Skew-symmetrical point forces. We assume now that the loading is given by two 
skew-symmetrical point forces again applied at a distance a behind crack tip and directed along the 
xa-axis, 

(P3>(2:i) = 0, lp:i}{xi)^-2F5{xi+a). (C.6) 

It is easy to show that all the formulae derived above still apply, however, they contain now the factor 
77, so that 



dxi IT y Xl Xl + a 



V^FJ — , (C.7) 
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rj — ( + ej arctann 4 / , — a < xi < 0, 

TT V a 

^i) - <( ^ (C.8) 

ri — fo + e) arctann . / , xi < — a, 

TT V ^1 

(cT23>(+^(.Ti)-77-./f (C.9) 



and the stress intensity factor is given by 



Km = m —F. (C.IO) 
V Tra 



Mode I and II. Symmetrical point forces. We consider now the vector 2D case and suppose 
that the loading is given as symmetrical point forces in both xi and X2 directions, namely 

{p{){xi)^-F^5{xi+a), M(2^i) = 0, (C.ll) 

{V2){xi) = -F26{xi+a), lp2\{xi) = Q. (C.12) 

For simplicity we assume that the Dundurs parameter d is equal to zero, so that the system of singular 
equations decouples 

45(^^^M:^ = fe)(xi), . = 1,2. (C.13) 
and it is possible to obtain results in closed form. Indeed, the inversion formula (jA.2|) leads to 



djujf-^ _ bFj i~r 1 

dxi TT y xi xi + 



(C.14) 



which, after integration, gives 



2bFj ^ , / xi 



arctanh w , — a < xi < 0, 



2bF, , I a 
^arctanh./ , xi < —a. 

TT \ Xl 

Finally, the tractions ahead of the crack tip are calculated from ([5T|) 



{a2,rHxi)^^.f^^. (C.16) 

TT \ Xl Xl + a 

The stress intensity factors are then obtained as 

Ki = J—F2, Ku = J—Fi. (C.17) 
V Tra \i TTa 

Mode I and II. Skew-symmetrical point forces. In this case the point forces, applied at a 
distance a behind the crack tip, have the same direction 

{pi){xi) = 0, lpiiixi) = -2FiS{xi+a), (C.18) 

{P2){xi) = 0, bsKa:^!) = -2F2S{xi + a). (C.19) 
Assuming again d = 0, the solution of the system of singular integral equations (|50p leads to 

dlmji-^ _ baFi , . . -bjF2Sixi+al (C.20) 



dxi TT Y xi Xl + a 
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aM(-) baF2 ^ ,uwK( ^ \ (^o^^ 

- ' ' h^Fib{xi + a), (<-^-21) 



dxi TT Y xi xi + a 

which after integration gives 



M(2^i) = 



TT \ a 



2baFi , / xi 

■ arctanh w , —a < xi < 0, 



TT y Xl 



2baFi , I a , ^ 

■arctanh./ h 07_f 2, cci < — a, 



(C.22) 



25aF2 , / 2;i 

arctanhw , — a < xi < 0, 

TT V a 

M(^i) 

I 

TT V ""^l 

The tractions ahead of the crack tip are obtained from ([5T|) and read 



■ arctanh . / 07^1, •■J^i < —cl. 



OTi 2 TT y + a 

(-..>^-H:.i) = - i5(^)bil(xi) . (C.24) 

axi 2 TT y Xl Xl + a 

The stress intensity factors are then obtained as 

Ki = a\—F2, Ku = a\ —Fi. (C.25) 
V Tra \ na 

The reader may be surprised by the couphng constant terms in the right-hand side of expressions (jC.22p , 
showing a discontinuity of the crack opening at the point of apphcation of the forces. However, these 
terms can be easily explained by making recourse to the Flamant solution for a half-plane loaded by 
concentrated forces Fi and F2 at its surface (Barber, 2002). 

Assuming a reference system centred at the point of application of the forces, the surface displacement 
are determined, up to an arbitrary constant, by 

, , Fi(l-j.)log|xi| F2(l-2i.)sign(xi) 
wii^i) = + , 

TT/i 4/i 

(C.26) 

, , Fi(l-2i.)sign(xi) Fi(l-z.)log|xi| 
U2[Xl) = , 

where fi and v are the shear modulus and the Poisson's ratio of the half-plane, respectively. 

We now apply this solution to the two half-planes making up the cracked bimaterial plane. Of 
course, we do not expect this solution to match exactly with the correct solution obtained above, since 
the boundary conditions for the crack problem differ from the Flamant problem along the interface 
joining the two materials. Nonetheless, let us for the time being make this analogy and calculate from 
the solution (jC.26|) what the crack opening would be. 

We consider only the coupling (diagonal) terms in the system (jC.26|) , which can explain the observed 
discontinuity in the crack opening. We also adjust the arbitrary constant in the Flamant solution in 
order to have the same displacement for the two crack faces at the crack tip. As a result, we obtain 

KKa^i) = F2 sign(-xi - a) (^^^ + ^^^) = W sign(-xi - a), (C.27) 

luij{xi) = -i^i sign(-xi - a) (^^^ + ^^^) = " W sign(-xi - a), (C.28) 
which is in perfect agreement with our solution (|C.22p . 
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D Identities involving convolution integrals 

Theorem D.l. Let >f{xi) be a function defined on R and possess the derivatives of order 1, ...,n. Assume 
also that the function (p{xi) is vanishing with all its derivatives as xi — > ±00. Then 

-00 J-oc~ y^{xi + (2:3 -6)^ -'-00 2:1 - 9a;r 

Proof. Consider the integral 

R roo ^ Qn 

^dCirf^, (D.2) 



-R 



and integrate by parts the inner integral, to obtain 

1 9"-V 



xi-^i a"-v 



/ J-00 [(:C,-6)2 + (x3-6)2]3/^ a^r'^^'^^- 

(D-3) 

The first term in (jP.Sp is zero and the second term can be rewritten changing the order of integration 
in the form 

- r i r Tf FITTT rw7^^^3 1 (xi - 6)^^rf6- (D.4) 

\J-R [{Xl - 6)2 + (2:3 - 6)2]-^/2 / Q^n 



Taking the limit as i? — > cx3. the inner integral in (jD.4p gives 

L [(.i^6)2+k-6)2F^'^^ = 

which concludes the proof. □ 

Theorem D.2. Let ip{xi ,x^) he a function defined on and possess the Fourier transform with respect 
to xi and X3, denoted by ^(/3, A). Then 

" r ^1-6 dy^,, , r r ^3 ~ 6 dip ^ 

.00 i-oo i^l - 6)2 + (X3 - 6)2 96 J -00 J -00 i^l 6)2 + (X3 - 6)2 96 (QJ) 

= 27r(p(x 1,2:3). 

Proof. Consider the Fourier transform A) and write 

^) = ^I^A^ • ^) + • A). (D.7) 

To conclude the proof, apply the inversion Fourier transform to both side of (jD.7p and use the first two 
formulae given in Tab. [T] □ 

Theorem D.3. Let ip{xi, x^) be a function defined on and possess the Fourier transform with respect 
X\ and x^, denoted by JpdS , X) . Then 

d6«6 = / / 7 ^ ^n , . rTI^"6a6- (D.8) 



00 •/ —00 



(xi - 6)2 + (2:3 - 6)2 96 7„oo 7-00 (2:1 - 6)2 + (X3 - 6)2 96 

Proof. Consider the Fourier transform Tp[f3, A) and write 

^ ,A<^(AA) = — ^•m/3,A). (D.9) 



To conclude the proof, apply the inversion Fourier transform to both side of (jD.7p and use the third and 
fourth formulae in Tab. [TJ □ 
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E Double Fourier transform inversion for the 3D case 

The following four integrals will be used in the analysis (Gradshteyn and Ryzhik, 1965): 



cos ax 



yW+x' 



zdx = Ko{al3), a> 0, Re/3 > 0, 



X sin ax , IT _„R 

— -dx = -e 

P^+x^ 2 



Xo(a/3)cos(6/3)d/3= - 



a > 0, Re^>0, 
1 



2 + P 



Re a > 0, 



e-"'^ cos{bf3)df3 



12 + 52 



, a > 0, 



(E.l) 

(E.2) 
(E.3) 
(E.4) 



where Ko{z) is the modified Bessel function of the second kind of zero order, satisfying the differential 
equation 

z^y" + zy' - z^y = 0. (E.5) 

In some sense, the last two integrals are the Fourier inversion formulae for the first two. 

We also shall make use of the fact that the Fourier transform of the convolution is given by the 
product of the Fourier transforms, which leads to 



J-^-^J-,-i[5(/3,A)-/(/3,A)] = 



00 pOO 



(E.6) 



— 00 J —OQ 



It is evident what to do with the factors /3 and A. Namely, these factors arise when taking the Fourier 
transform of the derivatives instead of the function itself: 



/3/(/3,A) =z^(/3,A), A/(/3,A) =*^(/3,A). 



(E.7) 



Analogous identities can be written for the factors 0^ , \^ and /3A involving higher order derivatives: 



^2/(/3,A) = -Q(/3,A), a2/(/3,A) = -Q(/3,A), /3A/(/3, A) 



dxidx^ 



(/3,A). (E., 



Taking into account the following identity: 



VF+A2 = 



/32 + A2 



we need to evaluate only the following three basic integrals: 



•^13 -^A 



-p-l-p-l 
-^13 -^A 



1 



A 

/32 + A2 
/32 + A2 



1 

47r2 
1 

1 

47J-2 



00 pOQ 



1 



00 pOO 



-00 v//?2 + A2 
A 



— 00 ^ — 00 
00 nOO 



^2 + A2 



ixiP —ixsX 



dpdX, 



00 ^ —00 



^2 + A2 



-?xi/3 — -ia^s A 



d(id\. 



Let us consider the first integral (jE.lOp . By (jE.ip we obtain 



-2x3 A 



dX 



cos(|a^3|A) 

yFTA2 - VFTA2 



dA = 2 ■ 



/^o(|x3|/3), Re/3>0, 
A'o(-|x3|/3), -Re/3>0. 



(E.9) 

(E.IO) 

(E.ll) 
(E.12) 

(E.13) 
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Making recourse to (jE.3[) , we conclude that 



^ / e-"^^Xo(|x3||/3|)d/3, 



Ko{\x3\l3)cos{xiP)dp, 



1 1 



(E.14) 



27r y^xl + xl 

We note the the function transforms into itself (multiplied by a constant factor) after double Fourier 
transform. 

Let us now consider the second integral (jE.lip . By means of (jE.2p we write 



/32 + A2 



-2isign(a;3) 



Asin(|a;3|A) 

;32 + A2 



d\ = — 7ri sign(a;3) 



-l^^l^, Re^>0, 



We can now use the formula (IE. 41) to obtain 



Re /3 > 0. 



(E.15) 



A 



/32 + A2 



-zsign(a;3) 
47r 

-isign(a;3) 



dl3, 



2tt 

-isign(a;3) |a;3| 



e-l^=^l'^cos(xi;3)d^, 



-l X3 



(E.16) 



2tt xf + x^ 2tt xf + x^ 

Again, the function transforms into itself (multiplied by a constant factor) after double Fourier transform. 
The last integral (|E.12p can be evaluated in a similar way and we have 



/3 



/32 + A2 



2tt x\ + xl 



(E.17) 



Equations (jE.14p , (|E.16p and (jE.17|) are the basic integrals we use to perform the double Fourier inversion 
for the 3D case in Sec. [5] The complete list of inversion formulae is given in Tab. [TJ 
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Function: Double Fourier inversion: 

/3 n- I xi dip 

I3(p 



13^ + 27rxl+xl dxi 

A I xs ^ dip 

\p 



/32 + A2 ^ 27r x\+xl 8x3 

P \ xi dip 

\ip ^' 



/32 + A2 ^ 27rx2+a;2 

A 1 i9(/j 



;32 + A2 27rx2+a:2 dxi 

1 „_ « 1 9<y9 

1 , I 1 

+ A2 ^ 27r ^a;2 + dxs 

1 „2- 1 1 



1 2- 1 1 <9V 



V^^+A^ 



A'^ 



v/^+A^ 
+ A2 • ^ 



27r 


V^l + 4 


1 


1 


27r 


y^xl+xl 


1 


1 


27r 


^/xl+xl 


1 


1 



l3Xip 



2^y^fT^ \dxl dxi 



Table 1: Double Fourier inversion formulae 
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